

/*
This do-file is the main analysis do-file for the zip level analysis
*/

* ++++++++++++++++++++++	
* ++++++++++++++++++++++	
*    SUM STATS TABLE
* ++++++++++++++++++++++
* ++++++++++++++++++++++	

* summary statistics table 
use "${data_derived}/zip_regression_data.dta", clear

* make data unique on zcta
keep if year == 2020 & week < 31
bys zcta : keep if _n ==1 

* keep only those zctas for which variates are non-missing
foreach var in men white black asian med_hh_inc total_mngmt_bus_sci_arts total_service_occ ///
	total_prod_trans total_age_below_18 total_age_18_24 total_age_25_34 /// 
	total_age_35_44 total_age_45_54 total_age_55_64 total_age_65_74 total_age_75_and_above ///
	hs_ged some_coll coll pop_density frac_broadband pop number_poi_total {
		keep if `var' !=. 
	}

* get means and SDs for high and low exposure zctas 
replace pop_density = pop_density * 1000000
gen var_name = ""
gen mean_low = . 
gen sd_low = . 
gen mean_high = . 
gen sd_high = . 
local i = 1
foreach var in men white black asian med_hh_inc total_mngmt_bus_sci_arts total_service_occ ///
	total_prod_trans total_age_below_18 total_age_18_24 total_age_25_34 /// 
	total_age_35_44 total_age_45_54 total_age_55_64 total_age_65_74 total_age_75_and_above ///
	hs_ged some_coll coll pop_density frac_broadband pop number_poi_total {
	
	quietly {
		su `var' [w=pop] if high_exp == 0, d
		replace var_name = "`var'" in `i'
		replace mean_low = `r(mean)' in `i'
		replace sd_low = `r(sd)' in `i'
		su `var' [w=pop] if high_exp == 1, d
		replace mean_high = `r(mean)' in `i'
		replace sd_high = `r(sd)' in `i'
		local ++i
	}
}

qui: count if high_exp == 0 
replace var_name = "n_zcta" in `i'
replace mean_low = `r(N)' in `i'
qui: count if high_exp == 1 
replace mean_high = `r(N)' in `i'

keep var_name mean_low sd_low mean_high sd_high
order var_name mean_low sd_low mean_high sd_high
keep if var_name != ""

export delimited "${tables}/zcta_sum_stats.csv", replace

* ++++++++++++++++++++++	
* ++++++++++++++++++++++	
*      REGRESSIONS
* ++++++++++++++++++++++
* ++++++++++++++++++++++	

* zip level regression data by week
use "${data_derived}/zip_regression_data.dta", clear
drop if week >=31

* interact high exposure indicator with week and year separately
* also interact high exposure indicator with week and year combined --> i.e. triple interact
qui: su week if year ==2020
local max_week = `r(max)'
foreach y in 2019 2020 {
	foreach var in high_exp {
		qui: gen `var'_`y' = `var' if year == `y'
		qui: replace `var'_`y' = 0 if year != `y'
		forval i=1/`max_week' {
			if `y' ==2020 {
				qui: gen `var'_`i' = `var' if week == `i'
				qui: replace `var'_`i' = 0 if week !=`i'
				qui: gen `var'_`i'_`y' = `var' if week == `i' & year == `y'
				qui: replace `var'_`i'_`y' = 0 if week !=`i' | year != `y'
				}
			}
		}		
	}	

* interact covariates with week FE
* this includes SCI weighted exposure to high-income places, high pop density places, etc.
qui: su week if year ==2020
forval i=1/`r(max)' {
	foreach var of varlist men white black asian med_hh_inc total_service_occ ///
		total_prod_trans total_mngmt_bus_sci_arts total_age_below_18 ///
		total_age_18_24 total_age_25_34 total_age_35_44 total_age_45_54 ///
		total_age_55_64 total_age_65_74 total_age_75_and_above hs_ged ///
		some_coll coll pop_density frac_broadband ///
		vent_nat_pc_inc vent_nat_pc_pop_dens vent_nat_pc_frc_urban ///
		vent_cty_pc_c_11_pc vent_cty_pc_c_13_pc vent_cty_pc_c_17_pc {
		qui: gen `var'_`i'_2020 = `var' if week == `i' & year ==2020 
		qui: replace `var'_`i'_2020 = 0 if week !=`i' | year !=2020
		}
	}	
	
* ++++++++++++++++++++++	
* social distancing
* ++++++++++++++++++++++

* regressions for fraction of devices home as outcome variable
order men_* white_* black_* asian_* med_hh_inc_* total_mngmt_bus_sci_arts_* ///
	total_service_occ_* total_prod_trans_* ///
	total_age_below_18* total_age_18_24* total_age_25_34* total_age_35_44* ///
	total_age_45_54* total_age_55_64* total_age_65_74* total_age_75_and_above* ///
	hs_ged_* some_coll_* ///
	coll_* pop_density_* frac_broadband_* ///
	vent_nat_pc_inc_* vent_nat_pc_pop_dens_* vent_nat_pc_frc_urban_*, last
order high_exp_*_2020, last
qui: su week if pct_devices_home !=. & year == 2020
local max_week_sd = `r(max)'
foreach var in pct_chg_dist_trvld pct_devices_home {
	reghdfe `var' high_exp_2-high_exp_`max_week_sd' high_exp high_exp_1 [w=pop], ///
		absorb(fips#week zcta) cluster(zcta)
		
	* save coefficient estimates	
	qui: cap gen x=.
	qui: gen c_`var'= . 
	qui: gen se_c_`var' =. 
	local u = 1
	forval i=1/`max_week_sd' {	
		qui: replace c_`var' = _b[high_exp_`i'] in `u'
		qui: replace se_c_`var' = _se[high_exp_`i'] in `u'
		qui: replace x = `u' in `u'
		local ++u
		}
	qui: gen ci_u_`var' = c_`var' + 1.96 * se_c_`var'
	qui: gen ci_l_`var' = c_`var' - 1.96 * se_c_`var'
}
		
		
* this time we allow for time-varying controls
qui: su week if pct_devices_home !=. & year == 2020
local max_week_sd = `r(max)'
foreach var in pct_chg_dist_trvld pct_devices_home {		
	reghdfe `var' high_exp_2-high_exp_`max_week_sd' high_exp high_exp_1 ///
		men_2_2020-men_`max_week_sd'_2020 men men_1_2020 ///
		white_2_2020-white_`max_week_sd'_2020 white white_1_2020 ///
		black_2_2020-black_`max_week_sd'_2020 black black_1_2020 ///
		asian_2_2020-asian_`max_week_sd'_2020 asian asian_1_2020 ///
		med_hh_inc_2_2020-med_hh_inc_`max_week_sd'_2020 med_hh_inc med_hh_inc_1_2020 ///
		total_mngmt_bus_sci_arts_2_2020-total_mngmt_bus_sci_arts_`max_week_sd'_2020 total_mngmt_bus_sci_arts total_mngmt_bus_sci_arts_1_2020 ///
		total_service_occ_2_2020-total_service_occ_`max_week_sd'_2020 total_service_occ total_service_occ_1_2020 ///
		total_prod_trans_2_2020-total_prod_trans_`max_week_sd'_2020 total_prod_trans total_prod_trans_1_2020 ///
		total_age_below_18_2_2020-total_age_below_18_`max_week_sd'_2020 total_age_below_18 total_age_below_18_1_2020 ///
		total_age_18_24_2_2020-total_age_18_24_`max_week_sd'_2020 total_age_18_24 total_age_18_24_1_2020 ///
		total_age_25_34_2_2020-total_age_25_34_`max_week_sd'_2020 total_age_25_34 total_age_25_34_1_2020 ///
		total_age_35_44_2_2020-total_age_35_44_`max_week_sd'_2020 total_age_35_44 total_age_35_44_1_2020 ///
		total_age_45_54_2_2020-total_age_45_54_`max_week_sd'_2020 total_age_45_54 total_age_45_54_1_2020 ///
		total_age_55_64_2_2020-total_age_55_64_`max_week_sd'_2020 total_age_55_64 total_age_55_64_1_2020 ///
		total_age_65_74_2_2020-total_age_65_74_`max_week_sd'_2020 total_age_65_74 total_age_65_74_1_2020 ///
		total_age_75_and_above_2_2020-total_age_75_and_above_`max_week_sd'_2020 total_age_75_and_above total_age_75_and_above_1_2020 ///
		hs_ged_2_2020-hs_ged_`max_week_sd'_2020 hs_ged hs_ged_1_2020 ///
		some_coll_2_2020-some_coll_`max_week_sd'_2020 some_coll some_coll_1_2020 ///
		coll_2_2020-coll_`max_week_sd'_2020 coll coll_1_2020 ///
		pop_density_2_2020-pop_density_`max_week_sd'_2020 pop_density pop_density_1_2020 ///
		frac_broadband_2_2020-frac_broadband_`max_week_sd'_2020 frac_broadband frac_broadband_1_2020 ///
		vent_nat_pc_inc_2_2020-vent_nat_pc_inc_`max_week_sd'_2020 vent_nat_pc_inc vent_nat_pc_inc_1_2020 ///
		vent_nat_pc_pop_dens_2_2020-vent_nat_pc_pop_dens_`max_week_sd'_2020 vent_nat_pc_pop_dens vent_nat_pc_pop_dens_1_2020 ///
		vent_nat_pc_frc_urban_2_2020-vent_nat_pc_frc_urban_`max_week_sd'_2020 vent_nat_pc_frc_urban vent_nat_pc_frc_urban_1_2020 ///
		[w=pop], absorb(fips#week zcta) cluster(zcta) tol(1e-5)
			
	* save coefficient estimates	
	qui: cap gen x=.
	qui: gen c_`var'_cov = . 
	qui: gen se_c_`var'_cov =. 
	local u = 1
	forval i=1/`max_week_sd' {	
		qui: replace c_`var'_cov = _b[high_exp_`i'] in `u'
		qui: replace se_c_`var'_cov = _se[high_exp_`i'] in `u'
		qui: replace x = `u' in `u'
		local ++u
		}
	qui: gen ci_u_`var'_cov = c_`var'_cov + 1.96 * se_c_`var'_cov
	qui: gen ci_l_`var'_cov = c_`var'_cov - 1.96 * se_c_`var'_cov		
	}

* ++++++++++++++++++++++
* POI
* ++++++++++++++++++++++	
	
* loop over regressions for different types of POIs
order men_* white_* black_* asian_* med_hh_inc_* total_mngmt_bus_sci_arts_* ///
	total_service_occ_* total_prod_trans_* ///
	total_age_below_18* total_age_18_24* total_age_25_34* total_age_35_44* ///
	total_age_45_54* total_age_55_64* total_age_65_74* total_age_75_and_above* ///
	hs_ged_* some_coll_* ///
	coll_* pop_density_* frac_broadband_* ///
	vent_nat_pc_inc_* vent_nat_pc_pop_dens_* vent_nat_pc_frc_urban_*, last
order high_exp_*_2020, last
qui: su week if l_num_visits_total !=. & year == 2020
local max_week_poi = `r(max)'
foreach var in total entertain fd_plces fd_stores health parks retail schools {		
	reghdfe l_num_visits_`var' high_exp_2-high_exp_`max_week_poi' high_exp high_exp_1 ///
		[w=number_poi_`var'], absorb(fips#week zcta) cluster(zcta)
		
	* save coefficient estimates	
	qui: cap gen x=.
	qui: gen c_l_num_visits_`var' = . 
	qui: gen se_c_l_num_visits_`var' =. 
	local u = 1
	forval i=1/`max_week_poi' {	
		qui: replace c_l_num_visits_`var' = _b[high_exp_`i'] in `u'
		qui: replace se_c_l_num_visits_`var' = _se[high_exp_`i'] in `u'
		qui: replace x = `u' in `u'
		local ++u
		}
	qui: gen ci_u_l_num_visits_`var' = c_l_num_visits_`var' + 1.96 * se_c_l_num_visits_`var'
	qui: gen ci_l_l_num_visits_`var' = c_l_num_visits_`var' - 1.96 * se_c_l_num_visits_`var'	
	}
	
* loop over regressions for different types of POIs
* this time we allow for time-varying controls
qui: su week if l_num_visits_total !=. & year == 2020
local max_week_poi = `r(max)'
foreach var in total entertain fd_plces fd_stores health parks retail schools {		
	reghdfe l_num_visits_`var' high_exp_2-high_exp_`max_week_poi' high_exp high_exp_1 ///
		men_2_2020-men_`max_week_poi'_2020 men men_1_2020 ///
		white_2_2020-white_`max_week_poi'_2020 white white_1_2020 ///
		black_2_2020-black_`max_week_poi'_2020 black black_1_2020 ///
		asian_2_2020-asian_`max_week_poi'_2020 asian asian_1_2020 ///
		med_hh_inc_2_2020-med_hh_inc_`max_week_poi'_2020 med_hh_inc med_hh_inc_1_2020 ///
		total_mngmt_bus_sci_arts_2_2020-total_mngmt_bus_sci_arts_`max_week_poi'_2020 ///
			total_mngmt_bus_sci_arts total_mngmt_bus_sci_arts_1_2020 ///
		total_service_occ_2_2020-total_service_occ_`max_week_poi'_2020 total_service_occ total_service_occ_1_2020 ///
		total_prod_trans_2_2020-total_prod_trans_`max_week_poi'_2020 total_prod_trans total_prod_trans_1_2020 ///
		total_age_below_18_2_2020-total_age_below_18_`max_week_poi'_2020 total_age_below_18 total_age_below_18_1_2020 ///
		total_age_18_24_2_2020-total_age_18_24_`max_week_poi'_2020 total_age_18_24 total_age_18_24_1_2020 ///
		total_age_25_34_2_2020-total_age_25_34_`max_week_poi'_2020 total_age_25_34 total_age_25_34_1_2020 ///
		total_age_35_44_2_2020-total_age_35_44_`max_week_poi'_2020 total_age_35_44 total_age_35_44_1_2020 ///
		total_age_45_54_2_2020-total_age_45_54_`max_week_poi'_2020 total_age_45_54 total_age_45_54_1_2020 ///
		total_age_55_64_2_2020-total_age_55_64_`max_week_poi'_2020 total_age_55_64 total_age_55_64_1_2020 ///
		total_age_65_74_2_2020-total_age_65_74_`max_week_poi'_2020 total_age_65_74 total_age_65_74_1_2020 ///
		total_age_75_and_above_2_2020-total_age_75_and_above_`max_week_poi'_2020 total_age_75_and_above total_age_75_and_above_1_2020 ///
		hs_ged_2_2020-hs_ged_`max_week_poi'_2020 hs_ged hs_ged_1_2020 ///
		some_coll_2_2020-some_coll_`max_week_poi'_2020 some_coll some_coll_1_2020 ///
		coll_2_2020-coll_`max_week_poi'_2020 coll coll_1_2020 ///
		pop_density_2_2020-pop_density_`max_week_poi'_2020 pop_density pop_density_1_2020 ///
		frac_broadband_2_2020-frac_broadband_`max_week_poi'_2020 frac_broadband frac_broadband_1_2020 ///
		vent_nat_pc_inc_2_2020-vent_nat_pc_inc_`max_week_poi'_2020 vent_nat_pc_inc vent_nat_pc_inc_1_2020 ///
		vent_nat_pc_pop_dens_2_2020-vent_nat_pc_pop_dens_`max_week_poi'_2020 vent_nat_pc_pop_dens vent_nat_pc_pop_dens_1_2020 ///
		vent_nat_pc_frc_urban_2_2020-vent_nat_pc_frc_urban_`max_week_poi'_2020 vent_nat_pc_frc_urban vent_nat_pc_frc_urban_1_2020 ///
		[w=number_poi_`var'], absorb(fips#week zcta) cluster(zcta) tol(1e-5)
		
	* save coefficient estimates	
	qui: cap gen x=.
	qui: gen c_l_num_visits_`var'_cov = . 
	qui: gen se_c_l_num_visits_`var'_cov =. 
	local u = 1
	forval i=1/`max_week_poi' {	
		qui: replace c_l_num_visits_`var'_cov = _b[high_exp_`i'] in `u'
		qui: replace se_c_l_num_visits_`var'_cov = _se[high_exp_`i'] in `u'
		qui: replace x = `u' in `u'
		local ++u
		}
	qui: gen ci_u_l_num_visits_`var'_cov = c_l_num_visits_`var'_cov + 1.96 * se_c_l_num_visits_`var'_cov
	qui: gen ci_l_l_num_visits_`var'_cov = c_l_num_visits_`var'_cov - 1.96 * se_c_l_num_visits_`var'_cov	
	}	
	
* ++++++++++++++++++++++	
* save the resulting data
* ++++++++++++++++++++++

keep c_* se_c_* ci_l_* ci_u_* x
egen miss = rowtotal(*)
keep if miss != 0 | _n ==1 

save "${data_derived}/zip_regression_results.dta", replace
